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Abstract. 

We present numerical simulations of axisymmetric magnetised massive tori around rotating black 
^-y holes taking into account the energy losses due to emission of neutrinos. A realistic equation of state 

is used which takes into account the energy losses due to dissociation of nuclei. The heating due 
to neutrino-antineutrino annihilation is not included. We study the cases of optically thick, semi- 
transparent, and optically thin to neutrino discs. We show that neutrino cooling does not change 
significantly the structure of accretion flow and the total energy release. The time scale of accretion 
is set by the torus angular momentum. Due to the lack of magnetic dynamo in our calculations, it is 
(— i ' the initial strength of magnetic field and its topology that determine the process of jet formation and 

its energetics. Extrapolation of our results gives the total energy released in the jet ~ 10 52 erg. This 
is sufficient to explain the hypernovae explosions associated with GRB 980425 and GRB 030329. 
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: INTRODUCTION 

-xt- : 

The long-soft Gamma Ray Bursts (GRB) are found predominantly in the blue light 
of galaxies 01 0, @], that is in the region where the most massive stars die. This is 
supported by the debated connection of GRB 980425 with SN 1998bw fl H @] and the 
detection of low-redshift (z = 0.1685; [7]) GRB 030329 and its associated supernovae, 
SN 2003dh tfsL&Ql]. Moreover, the afterglow spectra of many GRBs show broad SN 
spectral features. 

The most popular model of GRB central engine is based on the "failed supernova" 
scenario of stellar collapse, or "collapsar", where the iron core of progenitor forms a BH 
iflTll . If the progenitor is non-rotating then its collapse is likely to continue in a "silent" 
manner until the whole star is swallowed by the BH. If, however, the specific angular 
momentum in the equatorial part of the stellar envelope exceeds that of the last stable 
orbit of the BH then the collapse becomes highly anisotropic. While in the polar region 
it may proceed more or less uninhibited the equatorial layers form dense and massive 
accretion disk. The gravitational energy released in the disk can be very large, more then 
sufficient to stop the collapse of outer layers and drive GRB outflows, presumably in the 
polar direction where the mass density is much lower [12, H]. 



The last few years witnessed dramatic progress in general relativistic magnetohy- 
drodynamic (GRMHD) simulations of BH accretion systems. They revealed complex 
structure that can be decomposed into a disk, corona, disk wind, and highly magne- 
tised polar region that hosts a black hole jet driven by the Blandford-Znajek mechanism 
[18. IT4L 15, 16, [13]. However, all these simulations used very primitive equations of 



state (EOS) and neglected neutrino cooling. The first 2D GRMHD simulations in Kerr 
metrics with realistic EOS and neutrino cooling simulations have been carried out only 
very recently [19] but the physical time span of their computations was rather short, only 
~ 0.06 sec. In this paper we present the results of a similar study whose main aim is to 
investigate the effects of neutrino cooling processes on the disk structure and dynamics. 



PHYSICAL PROCESSES 



In these simulations we use realistic EOS [20], which also accounts for the photo- 
disintegration of nuclei. The cooling rates due to neutrino emission via annihilation of 
pairs, photo-production, and plasma mechanisms are taken from 112 ill , via URCA 
processes from B22D, and via synchrotron mechanism from j23\\ . 

We compare 3 different cases: 1) without cooling (model DL2A9NC); 2) cooling in 
optically thin regime (model DL2A9C); 3) cooling with incorporated opacity effects 
(model DL2A9MC). The first case is not very realistic in the GRB context and we 
consider it only as a reference model. The second case suits well tori with mass of few 
percent of the solar mass. It could be relevant for the models of neutron star-neutron star 
and BH-white dwarf mergers. The third case requires a torus of several solar masses. 

To describe the neutrino opacity we follow [20], namely we multiply the local loss 
rate by the factor e~ Tv , where T v = S v nl v is the neutrino optical depth, n = p/m p 
is the concentration of baryons, m p is the proton mass, and S v is the neutrino in- 
teraction cross-section 0, |25l]. Neglecting the pair loading at high densities S v = 



([K e + K p ]Y e + K n {\ —Y e ))Te v , where e v = 3.15kT, K e = OoA e /2m 2 , m e is the electron 
mass in MeV, Oq ~ 1.71 X 10 -44 cm 2 , A ; - = (cy +ca) 2 + \ {cy — ca) 2 , where c v and ca 
are the vector and axial-vector coupling constant for given neutrino species. A e ~ 2.2, 
K n = Ob(l + 3g 2 )/(16m 2 ) and K p = C7 o [4sin 4 w - 2sin 2 0M/ + (1 + 3£ 2 )/4]/(4m 2 ), 
where sin 2 Oyy — 0.231 and gA — —1.26 is the axial-vector coupling constant. The char- 
acteristic length-scale of neutrino absorption l v is difficult to calculate directly due to the 
complex and dynamic structure of solutions. Instead we use the following simplification 

l v = ^ ~ ^0.27 + 1 j r g , where r g = GMbh I c 2 is the gravitational radius of BH. 

At very high densities p > 8.1 x 10 g cm~ 3 the process of electron capture on 
nuclei becomes important. In order to take this into account we have utilised the results 
obtained in [26] and made a simple approximation which is 7% accurate. 



SIMULATION SETUP 

The simulations are carried out with an upwind conservative scheme that based on a 
linear Riemann solver and uses the constrained transport method to evolve magnetic 
field. The details of this numerical method and various tests are described in rf29l 
[3Qh . The neutrino cooling is introduced via the source term in the energy-momentum 
equation 



where u v is the four-velocity of plasma, g is the determinant of metric tensor, and 
S is the cooling rate as measured in the fluid frame. The gravitational attraction of 
BH is introduced via the Kerr metric in Kerr-Schild coordinates, {t,(j),r, 6}. The two- 
dimensional computational domain is (ro < r < r\) x (0 < 6 < it), where tq = (1 + 
Vl — a 2 /2) r g , there r g = 14.847 km and r\ = 200 r g = 2969 km. The total mass within 
the domain is small compared to the mass of BH (less then 25%) that allows us to ignore 
its self-gravity. The grid is uniform in 6 where it has 320 cells and almost uniform in 
log(r) where it has 459 cells, the linear cell size being the same in both directions. 

The initial solution describes an equilibrium torus with constant specific angular 
momentum Iq rf3~0. 321. The black hole mass Mbh = 10M Q and its rotation parameter 
a = 0.9. The torus mass M tor = 2.55M© and its specific angular momentum Iq = 2.8r g c = 
1.25 x 10 17 g cm 2 sec~ l . This setup corresponds to the collapse of a massive star when 
its core and inner envelope form a black hole and the outer envelope forms an accretion 
disk. The initial magnetic field is purely poloidal, thus the only non-zero component of 

its vector potential is A* °c W(r, 0) 3 , here W(r, 0) = A In iimzi!!£§ . The ratio of 

magnetic to gas pressure /3 = P m /P g does not exceed 3 x 10~ 2 . In the initial solution 
more then 80% of the total pressure comes from the radiation and only 20% from the 
degenerated matter. In these conditions the neutrino cooling is important. 



RESULTS 

After initial period of relaxation all our models show the same structure involving a 
thick disk, a highly magnetised jet and a wind from the disk. The fast neutrino cooling 
in the optically thin model DL2A9C leads to significant imbalance of the torus. The 
torus collapses to a new state where its thermal pressure becomes negligible. In the 
model DL2A9MC, where we take into account the opacity effects, the torus interior is 
optically thick to neutrinos. This model develops a slightly different structure where the 
hallo is squeezed as in model DL2A9MC but the central core is bigger and has a higher 
temperature (T ~ 2.2 x 10 11 K) compared to the model DL2A9C. Figure Q] shows the 
impact of neutrino cooling on the structure of the torus. Without cooling the torus keeps 
the same size and shape as in the initial solution. The strong cooling in model DL2A9C 
leads to the collapse of initial configuration with the peak density increasing 10-fold and 
reaching 10 13 g cm~ 3 . The moderate cooling in model DL2A9MC results in lower peak 
density, 3 x 10 12 g cm~ 3 . In DL2A9MC and DL2A9C the accretion flow is thinner than 
in DL2A9NC (see figCD). In DL2A9NC the highly magnetised jet region is much smaller 
and more variable than in DL2A9C and DL2A9MC. However, the global structure of the 
accretion flow does not change very much and the energy release rate in jets of all models 
is similar. In model DL2A9MC the efficiency, 

pT rT 

7] = / Etotdt/ / M BH c 2 dt, 
Jo Jo 

is about 0.003. In model DL2A9NC it is 3.8 times lower. 

At the start of simulations the central regions of models with cooling loose pressure 
support and begin to fall into the potential well. This leads to the oscillations similar to 



FIGURE 1. log 10 (p) (colour image) and magnetic field lines. The initial solution for all models is 
shown in the upper left panel. The other panels from left to right show solutions after 0.2075 sec; model 
DL2A9NC (second on the left), model DL2A9C (second on the right), model DL2A9MC (first on the 
right). 




FIGURE 2. Evolution of the mass-averaged torus radius. Dashed line - DL2A9NC, dot-dashed - 
DL2A9MC, solid - DL2A9C. 



those reported in [33]. However, we observe fast dumping of these oscillations (see fig. 
[2]) One reason is the strong radiative shocks driven into the outer layers of the torus which 
are slow to react to the loss of equilibrium. The other reason is Maxwell stresses. In the 
optically thin case (model DL2A9C) the initial amplitude of oscillations dr m / r™ = 0.05, 
where r m = J M rdm/ J M dm is the mass-averaged torus radius. The oscillations degrade 
quickly with the quality factor Q = \{drf + drf +l ) /2{dr? - drf +l )\ = 3.8, here drf 
is maximal deviation from equilibrium of i-th oscillation. The estimated relaxation 
time is t re [ ax = t osc x(3~ 0.027 sec. In model DL2A9MC with less effective cooling 
the initial amplitude of oscillations is higher, dr m /r m = 0.15, but the quality factor is 
lower, Q = 1.6. After few oscillations the amplitude decreases to dr" 1 /r m = 0.05, the 
effective viscosity decreases and the quality factor grows to Q = 6.7. The corresponding 
relaxation time t re \ ax ~ 0.063 sec. When neutrino cooling is not included the oscillations 
do not appear at all and the mass-averaged radius of the torus increases in time due to 
the development of wind. 



TABLE 1. Main results. 



model 


I 


a 


M BH 




F tot 


n 


DL2A9NC 


2.8 


0.9 


1.3937 


0.0189 


1.8530 


0.00074 


DL2A9C 


2.8 


0.9 


0.8588 


0.0134 


2.9039 


0.0019 


DL2A9MC 


2.8 


0.9 


0.512 


0.0276 


2.577 


0.0028 



DISCUSSION AND CONCLUSIONS 

In all our models we observe the development of a Poynting-dominated jet from BH. 
The neutrino cooling processes do not have a significant effect on the jet power. In fact 
it appears to be less disrupted due to interactions with disk wind and corona. On the 
other hand, the disk wind weakens when the cooling is included. The estimated total 
energy release during the whole period of accretion E tot ps r\M torus c 2 ~ 10 52 erg. This 
can explain the energetics of hypernovae explosions associated with GRB 980425 and 
GRB 030329. The studies of magnetic acceleration of Poynting-dominated jets show 
that high terminal Lorentz factors can be achieved together with effective conversion of 
electromagnetic energy into the bulk motion energy of the jet rf34l [351. Thus, this is a 
promising model of central engines of Gamma Ray Bursts. 

Our model is complimentary to the model of magnetar driven GRBs and hypernovae, 
e.g. ll3~6ll . Both these models explain the origin of collimated GRB jets and power 
of hypernovae. Further investigation is required to distinguish these two scenario and 
determine their specific observational signatures. At the moment we can only say that 
the chemical composition of the jet will different. In the collapsar case we expect highly 
Poynting dominated e + e~ jets from onset whereas in the magnetar model the jets are 
initially strongly contaminated by baryons. 

The disc oscillations observed in our simulations are the product of initial setup be- 
cause the effects of neutrino cooling have not been included in the initial equilibrium 
solution. They are quickly dumped and thus unlikely to be excited in realistic astrophys- 
ical conditions. This conclusion is supported by the simulations where the formation of 
accretion disks is followed from the onset of stellar collapse fl2L [l31. 
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